using System;
using System.Collections.Generic;
using System.Text;
using Algorithms.Numerical.Impl;

namespace Algorithms.Numerical
{
    public class MatrixFactory
    {
        public static MatrixFactory Instance
        {
            get
            {
                if (instance == null)
                    instance = new MatrixFactory();
                return instance;
            }
        }
            
        public Matrix MakeZero(int height, int width)
        {
            return new MatrixImpl(height, width);
        }

        public Matrix MakeIdentity(int height, int width)
        {
            Matrix m = this.MakeZero(height, width);
            for (int i = 0; i < Math.Min(height, width); i++)
                m[i, i] = 1;
            return m;
        }

        public Matrix Make(double[,] m)
        {
            return new MatrixImpl(m);
        }

        /*
         * Return matrix column conatenation:
         *  R = A|B
         */
        public Matrix ColumnsJoin(Matrix a, Matrix b)
        {
            if (a.Height != b.Height)
                throw new ArgumentException();
            Matrix r = MakeZero(a.Height, a.Width + b.Width);
            r.SubmatrixView(0, 0, a.Height, a.Width).Assign(a);
            r.SubmatrixView(0, a.Width, a.Height, b.Width).Assign(b);
            return r;
        }

        /*
         * Return matrix row conatenation:
         *      A
         *  R = -
         *      B
         */
        public Matrix RowsJoin(Matrix a, Matrix b)
        {
            if (a.Width != b.Width)
                throw new ArgumentException();
            Matrix r = MakeZero(a.Height + b.Height, a.Width);
            r.SubmatrixView(0, 0, a.Height, a.Width).Assign(a);
            r.SubmatrixView(a.Height, 0, b.Height, a.Width).Assign(b);
            return r;
        }

        /*
         * Compose matrix of matrices. Matrices in mm must have proper dimensions.
         */
        public Matrix Compose(Matrix[,] mm)
        {
            // TODO Speed could be improved. Now a lot of intermediate objects is created.

            Matrix[] a = new Matrix[mm.GetLength(0)];
            for (int y = 0; y < mm.GetLength(0); y++)
            {
                a[y] = mm[y, 0];
                for (int x = 1; x < mm.GetLength(1); x++)
                    a[y] = ColumnsJoin(a[y], mm[y, x]);
            }
            Matrix r = a[0];
            for (int y = 1; y < a.Length; y++)
                r = RowsJoin(r, a[y]);
            return r;
        }

        public Matrix Mult(double scale, Matrix a, bool assignToA)
        {
            Matrix r = a;
            if (!assignToA)
                r = MakeZero(a.Height, a.Width);
            for (int y = a.Height; --y >= 0 ; )
                for (int x = a.Width; --x >= 0; )
                    r[y, x] = scale * a[y, x];
            return r;
        }

        /*
         * R = alpha * a * b
         */
        public Matrix Mult(double alpha, Matrix a, Matrix b)
        {
            return Mult(alpha, a, b, 0, null, false );
        }

        /*
         * R = alpha * a * b + beta * c
         * 
         * c can be null. In such case zero matrix is used.
         */
        public Matrix Mult( double alpha, Matrix a, Matrix b, double beta, Matrix c, bool assignToC)
        {
            if (c != null && (a.Height != c.Height || b.Width != c.Width))
                throw new ArgumentException();
            if (a.Width != b.Height)
                throw new ArgumentException();
            if (c == null)
            {
                c = MakeZero(a.Height, b.Width);
                assignToC = true;
            }
            Matrix r = c;
            if ( ! assignToC )
                r = MakeZero(a.Height, b.Width);

            for( int y = a.Height; --y >= 0; )
                for (int x = b.Width; --x >= 0; )
                {
                    double s = 0;
                    for (int k = a.Width; --k >= 0; )
                        s += alpha * a[y, k] * b[k, x];
                    r[y, x] = s + beta * c[y, x];
                }
            return r;
        }

        public Matrix Sum(Matrix a, Matrix b, bool assignToA)
        {
            if (a.Width != b.Width || a.Height != b.Height)
                throw new ArgumentException();
            Matrix r = a;
            if ( ! assignToA )
                r = MakeZero(a.Height, a.Width);
            for (int y = a.Height; --y >= 0; )
                for (int x = b.Width; --x >= 0; )
                    r[y, x] = a[y, x] + b[y, x];
            return r;
        }

        
        






        private MatrixFactory()
        {
        }
        private static MatrixFactory instance = null;
    }
}
